home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / zhpr2.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  8.4 KB  |  255 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE ZHPR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, AP )
  5. *     .. Scalar Arguments ..
  6.       COMPLEX*16         ALPHA
  7.       INTEGER            INCX, INCY, N
  8.       CHARACTER*1        UPLO
  9. *     .. Array Arguments ..
  10.       COMPLEX*16         AP( * ), X( * ), Y( * )
  11. *     ..
  12. *
  13. *  Purpose
  14. *  =======
  15. *
  16. *  ZHPR2  performs the hermitian rank 2 operation
  17. *
  18. *     A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
  19. *
  20. *  where alpha is a scalar, x and y are n element vectors and A is an
  21. *  n by n hermitian matrix, supplied in packed form.
  22. *
  23. *  Parameters
  24. *  ==========
  25. *
  26. *  UPLO   - CHARACTER*1.
  27. *           On entry, UPLO specifies whether the upper or lower
  28. *           triangular part of the matrix A is supplied in the packed
  29. *           array AP as follows:
  30. *
  31. *              UPLO = 'U' or 'u'   The upper triangular part of A is
  32. *                                  supplied in AP.
  33. *
  34. *              UPLO = 'L' or 'l'   The lower triangular part of A is
  35. *                                  supplied in AP.
  36. *
  37. *           Unchanged on exit.
  38. *
  39. *  N      - INTEGER.
  40. *           On entry, N specifies the order of the matrix A.
  41. *           N must be at least zero.
  42. *           Unchanged on exit.
  43. *
  44. *  ALPHA  - COMPLEX*16      .
  45. *           On entry, ALPHA specifies the scalar alpha.
  46. *           Unchanged on exit.
  47. *
  48. *  X      - COMPLEX*16       array of dimension at least
  49. *           ( 1 + ( n - 1 )*abs( INCX ) ).
  50. *           Before entry, the incremented array X must contain the n
  51. *           element vector x.
  52. *           Unchanged on exit.
  53. *
  54. *  INCX   - INTEGER.
  55. *           On entry, INCX specifies the increment for the elements of
  56. *           X. INCX must not be zero.
  57. *           Unchanged on exit.
  58. *
  59. *  Y      - COMPLEX*16       array of dimension at least
  60. *           ( 1 + ( n - 1 )*abs( INCY ) ).
  61. *           Before entry, the incremented array Y must contain the n
  62. *           element vector y.
  63. *           Unchanged on exit.
  64. *
  65. *  INCY   - INTEGER.
  66. *           On entry, INCY specifies the increment for the elements of
  67. *           Y. INCY must not be zero.
  68. *           Unchanged on exit.
  69. *
  70. *  AP     - COMPLEX*16       array of DIMENSION at least
  71. *           ( ( n*( n + 1 ) )/2 ).
  72. *           Before entry with  UPLO = 'U' or 'u', the array AP must
  73. *           contain the upper triangular part of the hermitian matrix
  74. *           packed sequentially, column by column, so that AP( 1 )
  75. *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
  76. *           and a( 2, 2 ) respectively, and so on. On exit, the array
  77. *           AP is overwritten by the upper triangular part of the
  78. *           updated matrix.
  79. *           Before entry with UPLO = 'L' or 'l', the array AP must
  80. *           contain the lower triangular part of the hermitian matrix
  81. *           packed sequentially, column by column, so that AP( 1 )
  82. *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
  83. *           and a( 3, 1 ) respectively, and so on. On exit, the array
  84. *           AP is overwritten by the lower triangular part of the
  85. *           updated matrix.
  86. *           Note that the imaginary parts of the diagonal elements need
  87. *           not be set, they are assumed to be zero, and on exit they
  88. *           are set to zero.
  89. *
  90. *
  91. *  Level 2 Blas routine.
  92. *
  93. *  -- Written on 22-October-1986.
  94. *     Jack Dongarra, Argonne National Lab.
  95. *     Jeremy Du Croz, Nag Central Office.
  96. *     Sven Hammarling, Nag Central Office.
  97. *     Richard Hanson, Sandia National Labs.
  98. *
  99. *
  100. *     .. Parameters ..
  101.       COMPLEX*16         ZERO
  102.       PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
  103. *     .. Local Scalars ..
  104.       COMPLEX*16         TEMP1, TEMP2
  105.       INTEGER            I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY
  106. *     .. External Functions ..
  107.       LOGICAL            LSAME
  108.       EXTERNAL           LSAME
  109. *     .. External Subroutines ..
  110.       EXTERNAL           XERBLA
  111. *     .. Intrinsic Functions ..
  112.       INTRINSIC          DCONJG, DBLE
  113. *     ..
  114. *     .. Executable Statements ..
  115. *
  116. *     Test the input parameters.
  117. *
  118.       INFO = 0
  119.       IF     ( .NOT.LSAME( UPLO, 'U' ).AND.
  120.      $         .NOT.LSAME( UPLO, 'L' )      )THEN
  121.          INFO = 1
  122.       ELSE IF( N.LT.0 )THEN
  123.          INFO = 2
  124.       ELSE IF( INCX.EQ.0 )THEN
  125.          INFO = 5
  126.       ELSE IF( INCY.EQ.0 )THEN
  127.          INFO = 7
  128.       END IF
  129.       IF( INFO.NE.0 )THEN
  130.          CALL XERBLA( 'ZHPR2 ', INFO )
  131.          RETURN
  132.       END IF
  133. *
  134. *     Quick return if possible.
  135. *
  136.       IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
  137.      $   RETURN
  138. *
  139. *     Set up the start points in X and Y if the increments are not both
  140. *     unity.
  141. *
  142.       IF( ( INCX.NE.1 ).OR.( INCY.NE.1 ) )THEN
  143.          IF( INCX.GT.0 )THEN
  144.             KX = 1
  145.          ELSE
  146.             KX = 1 - ( N - 1 )*INCX
  147.          END IF
  148.          IF( INCY.GT.0 )THEN
  149.             KY = 1
  150.          ELSE
  151.             KY = 1 - ( N - 1 )*INCY
  152.          END IF
  153.          JX = KX
  154.          JY = KY
  155.       END IF
  156. *
  157. *     Start the operations. In this version the elements of the array AP
  158. *     are accessed sequentially with one pass through AP.
  159. *
  160.       KK = 1
  161.       IF( LSAME( UPLO, 'U' ) )THEN
  162. *
  163. *        Form  A  when upper triangle is stored in AP.
  164. *
  165.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  166.             DO 20, J = 1, N
  167.                IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN
  168.                   TEMP1 = ALPHA*DCONJG( Y( J ) )
  169.                   TEMP2 = DCONJG( ALPHA*X( J ) )
  170.                   K     = KK
  171.                   DO 10, I = 1, J - 1
  172.                      AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2
  173.                      K       = K       + 1
  174.    10             CONTINUE
  175.                   AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) ) +
  176.      $                               DBLE( X( J )*TEMP1 + Y( J )*TEMP2 )
  177.                ELSE
  178.                   AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) )
  179.                END IF
  180.                KK = KK + J
  181.    20       CONTINUE
  182.          ELSE
  183.             DO 40, J = 1, N
  184.                IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
  185.                   TEMP1 = ALPHA*DCONJG( Y( JY ) )
  186.                   TEMP2 = DCONJG( ALPHA*X( JX ) )
  187.                   IX    = KX
  188.                   IY    = KY
  189.                   DO 30, K = KK, KK + J - 2
  190.                      AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2
  191.                      IX      = IX      + INCX
  192.                      IY      = IY      + INCY
  193.    30             CONTINUE
  194.                   AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) ) +
  195.      $                               DBLE( X( JX )*TEMP1 +
  196.      $                                     Y( JY )*TEMP2 )
  197.                ELSE
  198.                   AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) )
  199.                END IF
  200.                JX = JX + INCX
  201.                JY = JY + INCY
  202.                KK = KK + J
  203.    40       CONTINUE
  204.          END IF
  205.       ELSE
  206. *
  207. *        Form  A  when lower triangle is stored in AP.
  208. *
  209.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  210.             DO 60, J = 1, N
  211.                IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN
  212.                   TEMP1   = ALPHA*DCONJG( Y( J ) )
  213.                   TEMP2   = DCONJG( ALPHA*X( J ) )
  214.                   AP( KK ) = DBLE( AP( KK ) ) +
  215.      $                       DBLE( X( J )*TEMP1 + Y( J )*TEMP2 )
  216.                   K        = KK               + 1
  217.                   DO 50, I = J + 1, N
  218.                      AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2
  219.                      K       = K       + 1
  220.    50             CONTINUE
  221.                ELSE
  222.                   AP( KK ) = DBLE( AP( KK ) )
  223.                END IF
  224.                KK = KK + N - J + 1
  225.    60       CONTINUE
  226.          ELSE
  227.             DO 80, J = 1, N
  228.                IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
  229.                   TEMP1    = ALPHA*DCONJG( Y( JY ) )
  230.                   TEMP2    = DCONJG( ALPHA*X( JX ) )
  231.                   AP( KK ) = DBLE( AP( KK ) ) +
  232.      $                       DBLE( X( JX )*TEMP1 + Y( JY )*TEMP2 )
  233.                   IX       = JX
  234.                   IY       = JY
  235.                   DO 70, K = KK + 1, KK + N - J
  236.                      IX      = IX      + INCX
  237.                      IY      = IY      + INCY
  238.                      AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2
  239.    70             CONTINUE
  240.                ELSE
  241.                   AP( KK ) = DBLE( AP( KK ) )
  242.                END IF
  243.                JX = JX + INCX
  244.                JY = JY + INCY
  245.                KK = KK + N - J + 1
  246.    80       CONTINUE
  247.          END IF
  248.       END IF
  249. *
  250.       RETURN
  251. *
  252. *     End of ZHPR2 .
  253. *
  254.       END
  255.